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ABSTRACT 


Computer Simulation techniques were used to determine 
equilibrium positions and binding energies of inert gas atoms 
implanted in a tungsten crystal and to investigate the po- 
tential wells around these equilibrium positions in both per- 
fect lattices and relaxed lattices. Stable positions were 
found for inert gas interstitials near lattice atoms in the 
Paine aid tourth layers of the crystal. Interstitial posi- 
tions near atoms in the first and second layers of the crys- 
tal appeared to be unstable if they exist at all. Asa 
result of potential well studies, it was concluded that the 
mechanism associated with equilibrium position formation was 
a combination of local liquefaction of the lattice structure 
and interaction of the interstitial with lattice atoms. 
Equilibrium positions were found to be ill-defined regions 
in the general (110) direction. The binding energy deter- 
mined for an interstitial site near a lattice atom in the 
third layer of the crystal was in excellent agreement with 


experimental results. 
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I. INTRODUCTION 


In 1968 Kornelsen and Sinha [1] of the National Research 
Council of Canada published results of radiation-damage ex- 
periments performed on tungsten. In these experiments a 
tungsten surface was bombarded with ions of neon, argon, 
krypton, and xenon respectively. Then, while the tungsten 
was heated, gas desorption rates were measured as the gas 
evolved from the crystal. The resultant desorption spectrum 
was interpreted to yield a binding energy spectrum for the 
trapped particles. In 1970 Professor Don E. Harrison, Jr., 
of the Naval Postgraduate School undertook the modeling of 
these experiments utilizing computer simulation techniques 
in order to provide a means for interpretation of their re- 
sults. Two successive thesis research efforts [2,3] have 
been specifically directed toward the investigation of inert 
gas implantation in a tungsten crystal. It was anticipated 
that a corollary to the successful computer simulation of 
this problem might possibly be an improved understanding of 
the interstitial Stor stabilization mechanism in tungsten, 
with more general application to other materials. 

The investigations reported in this paper are a con- 
tinuation of work begun by Vine [2] and Tankovich [3] under 
Harrison's supervision. The simulation procedures followed 
were a combination of static and dynamic MEE. to the 


problem. The static portion of the problem entailed 





interstitial implantation of an inert gas atom in a tung- 
sten crystal and the subsequent relaxation of the crystal 
until the equilibrium position of the interstitial could be 
ascertained. In the dynamic portion of the problem de- 
creasing amounts of energy were imparted to the interstitial 
in its equilibrium.position until the minimum energy, and 
direction, which still allowed the interstitial to escape 
from the crystal was determined (1.e., the binding energy 

of the particle). The binding energies thus determined 
provided a basis for comparison with the results of Kornelsen 


and Sinha. 





Peden) OF THE PROBLEM 


A. THE LATTICE DYNAMICS PROBLEM IN THE COMPUTER 

For a little over a decade high speed digital computers 
have exhibited their usefulness as tools for investigation 
of physical systems. Specifically, the inherent periodicity 
and order of crystal lattices have made the study of lattice 
dynamics particularly adaptable to Me EIE by computer 
ui Ton "Lt is a relatively simple matter to “construct” 
a crystal of the desired body-centered cubic or face-centered 
cubic structure in the computer. Various types of point de- 
fects can also be "created" in the crystal with relative 
ease. A vacancy is obtained by simply removing an atom from 
the crystal, while interstitials can be created by implanting 
an additional atom in the crystal. Two types of interstitials 
have been used in investigations of lattice dynamics, atoms 
identical to the crystal lattice atoms ("self-interstitials") 
and atoms different from the crystal lattice atoms 
(interstitials). 

Although a crystal lattice containing a point defect can 
be easily represented in a computer, modeling of the dynamics, 
which allows alterations of the crystal structure resulting 
from the presence of the point defect, is more involved. 

The dynamic portion of the problem of computer modeling of 
lattices can be characterized by four key decisions which 


must be made. 





Perse, 23t.is necessary to find some mathematical re- 
lationship to govern the interaction among the atoms of the 
crystal. This interaction is a complex many body problem 
which is nearly always approximated in computer simulations 
as a sum of appropriate two body interactions. With this 
approximation it is then possible to represent these two 
body interactions by some type of potential function, which, 
in turn, can be used to determine forces on individual atoms. 

Secondly, since the number of calculations required is 
directly related to the number of atoms in the crystal, a 
judicious choice of crystal size must be made. Large crystals 
would give more accurate results, but would also require more 
computational time. 

A third key problem which must be solved results from 
Dae smeabselity of a digital computer to perform a direct in- 
tegration. Since the lattice dynamics problem is most often 
directed at a determination of lattice atom positions after 
some type of interaction, an integration of the equations of 
motion is required. A choice of the numerical integration 
technique to be used must therefore be made. 

Last, and probably most important, since an iterative 
process is used to integrate the equations of motion, the 
length of the time interval of the iteration must be chosen 
so that the force variations within this time interval are 
small and, consequently, stability of the system is insured. 
At first glance, this suggests that only very small timesteps 


should be taken, but this would require excessive computer 
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time. Consequently, some variational method of timestep 
determination would be optimum. 
B. HISTORICAL DEVELOPMENT OF COMPUTER SIMULATION OF 

LATTICE DYNAMICS 

The four basic problems of computer simulation of lattice 
dynamics have been solved in various ways in previous 
investigations. 

In 1960 Gibson, Goland, Milgram, and Vineyard [4] 
(referred to hereafter as Gibson, et.al) of the Brookhaven 
National Laboratory published the first complete statement 
of computer simulation procedure as a tool for investigation 
of radiation damage of crystal materials. Due to the com- 
plexity of the radiation damage problem and the inadequacy 
of analytical methods as a means of analysis of the damage 
processes, Gibson, et.al. turned to a numerical integration 
of the problem utilizing high speed computers. Their pio- 
neering work gives insight into such pertinent subjects as 
crystal size, choice of potential functions, computational 
methods for solving the equations of motion, and timestep 
determination. 

Specifically, their research utilized copper as the 
crystal material since its relatively simple face-centered 
cubic structure and its widespread use in actual experimen- 
tation made it particularly appropriate for initial investi- 
gation. A Born-Mayer repulsive potential was used to describe 
atom-atom interaction, and a constant cohesive force was ap- 


plied to each atom on the crystal surface to balance the 
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Born-Mayer repulsion. A central difference iterative 
procedure was used to integrate the equations of motion. 

The procedure used for timestep determination is of parti- 
cular significance since the principles involved are still 
the governing criteria for choice of timestep duration. The 
fundamental observation was made that the greatest stress 

on the crystal system was a result of the strongest atom- 
atom interaction in the system. This interaction was there- 
fore chosen as the basis for timestep determination. Simply, 
the timestep duration was chosen to be inversly proportional 
to the velocity created by the strongest atom-atom 
interaction. 

In addition to verifying the applicability of computer 
Simulation techniques to radiation damage studies and pro- 
viding specific information on collision chains and focusing 
phenomena in crystallites, the work of Gibson, et.al. [4] 
determined that the only stable configuration for self- 
interstitials in a face-centered cubic crystal was the (100? 
split configuration. The split configuration implies that 
the interstitial causes a lattice atom to split away from 
its normal lattice position, and then both atoms, the inter- 
stitial and the lattice atom, share (or split) the normal 
lattice site when equilibrium is reached. 

Johnson and Brown [5] confirmed the split configuration 
as the only stable interstitial position in their studies 
of copper utilizing a Born-Mayer repulsive force between 
nearest neighbors and an elastic continuum containing the 


remainder of the atoms of the crystal.  Erginsoy, Vineyard, 
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and Englert [6] (referred to hereafter as Erginsoy, et.al.) 
extended these calculations to the body-centered cubic case 
by performing investigations using « iron. e ст 
techniques paralleled those of the Brookhaven group except 
for the choice of a potential function. After experimen- 
tation with a Born-Mayer potential and a Morse potential 
with parameters derived by Girifalco and Weizer [7], 
Erginsoy, et.al.[6] settled on a combination potential which 
utilized an exponentially screened Coulomb potential for 
close approaches, a Born-Mayer potential for first nearest 
neighbor, and a Morse or modified Morse potential for 
higher order neighbors. This research verified the split 
configuration as the only stable configuration for body- 
centered cubic structures, but the orientation of the split 
was found to be in a (110) plane, (vice a (100) plane as in 
the face-centered cubic case). R.A. Johnson [8] confirmed 
the split interstitial orientation in similar work with « 
iron, Vanadium and tungsten. 

In their work on collisions between a copper atom and 
a copper lattice, Gay and Harrison [9] introduced the average 
force CNET integrating the equations of motion. A 
complete description of this procedure has been reported by 
Gay, Effron, and Harrison [10]. In two more recent efforts 
Johnson and Wilson [11,12] determined new potential functions 
for use in face-centered cubic and body-centered cubic defect 
calculations and published results of defect calculations of 


heme various metals. 
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Parallel to the computer simulation efforts described 
above, Kornelsen [13,14] and Kornelsen and Sinha [1,15] 
have published the results of extensive experimentation 
involving the interaction of inert gas ions with tungsten. 
Under Harrison's Supervision, Vine [2] initiated attempts 
to devise a computer model of Kornelsen's experiments on 
neon implantation in tungsten. He attempted to develop a 
relationship between equilibrium positions of tungsten 
interstitials in a tungsten crystal determined, first, by 
using a Born-Mayer repulsive potential to describe tungsten 
interstitial interaction with a tungsten lattice then, by 
using the same tungsten-tungsten composite potential which 
was assumed between atomsof the tungsten crystal. This re- 
lationship would then have been applied to neon equilibrium 
positions derived from a repulsive potential to provide a 
comparison with Kornelsen's data. These efforts met with 
little success. 

Follow-on work by Tankovich [3] utilized the static/ 
dynamic approach described in more detail in Section III-B. 
The potential function used in these investigations was a 
composite Born-Mayer and Morse potential joined by the best 
cubic fit in the area of intersection which had previously 
been developed by Harrison and Moore [16]. Static program 
runs confirmed the (110) split interstitial for heliun, 
neon, argon, Krypton, and xenon point defects in a tungsten 
crystal for possible interstitial sites in planes three 


through six of the ten plane crystal used. Preliminary 
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Сиера е testing of an argon equilibrium position in plane 
four of the crystal was begun with limited success. 

Of particular significance in Tankovich's work [3] was 
the introduction of a timestep decrementing process into the 
static program. In previous investigations the timestep 
duration was chosen at the outset of the problem and remained 
constante through all computations. At this author's sug- 
gestion, a timestep decrementing process was devised which 
allowed a more rapid approach to the equilibrium positions 
with a concomitant saving of computer time required for 
computation. (See Section III-A-3 for a discussion of the 


timestep.) 
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III. THE SIMULATION PROCEDURE 


A. THE MODEL 
I Ihe crystal 

Two factors are of primary importance in determining 
the size of the crystal to be used in computer simulation 
investigations. First, the crystal must be large enough to 
provide realistic results, and, second, the crystal must be 
as small as possible in order to minimize computer time re- 
dared for calculation of atom-atom interactions. 

After experimenting with crystals of various sizes, 
Tankovich [3] determined that a crystal with ten planes of 
atoms in each coordinate direction, which is equivalent to 
five unit cells, (10 x 10 x 10) was suitable for static in- 
vestigations with tungsten. The same crystal was consequently 
used in the investigations reported in this paper. The lat- 
tice unit, or distance between adjacent (100) planes of 
atoms, for a body-centered cubic tungsten crystal is 1.582. 
All distances in the computer simulation were measured in 
lattice units. The lattice constant for the tungsten crystal 
is 3.162 (or two lattice units), and the nearest neighbor 
За засели Tattieesunits. 

The same numbering sequence for atoms of the crystal 
that was employed by Tankovich [3] was used in these investi- 
gation. (See Figure 1.) Atom number one was always assigned 
to the interstitial atom. A rectangular coordinate system 


was placed on the upper, left hand, front face of the crystal. 
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The positive "x" direction was chosen to the right from the 
origin, the positive "y" direction was down, and the positive 
"2" direction was to the rear. Atoms were numbered consecu- 
tively, beginning with atom number 2 at the origin and con- 
tinuing until all atoms of the Y = O plane had been numbered. 
This same numbering sequence was then followed for each Y 
plane of the crystal until all 250 atoms of the crystal had 
been numbered. Figure 1 shows the numbering of the atoms 
in the Y = 0 and Y = 1 planes. The numbering of atoms in the 
remainder of the planes follows the same pattern. 

For the dynamic program it was felt that a smaller 
crystal could be used and still yield meaningful results. 
The rationale for this determination was based on the go-no 
go (i.e., escape or no escape) character of the dynamic pro- 
gram. The dynamic program essentially provides an answer to 
this question - Will an ee el atom escape from the 
crystal if it is given a specific energy and directed in a 
specific direction? If the energy dissipation mechanism 
provided by collisions with, and close approaches to, the 
lattice atoms along the path traveled is great enough, the 
interstitial will remain in the crystal. If these energy 
dissipation mechanisms are not large enough to overcome the 
kinetic energy of the interstitial, the interstitial will 
escape. Since the atoms of the crystal along this line of 
motion must provide the mechanism to dissipate the inter- 
stitial's kinetic energy, for escape to be prevented, only 
atoms in the horizontal planes above the interstitial and 


atoms in vertical planes within a few lattice units of the 
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interstitial will affect the possible escape. The remainder 
of the atoms of the crystal would not have time to react 
with the interstitial or affect lu movement. As a minimum, 
the smaller crystal could be used to eliminate excessively 
high or low energies from consideration at a considerable 
savings of computer time and, instead, give a limited range 
of energies to be checked by dynamic runs using the entire 
crystal. 

To implement this procedure SUBROUTINE PLUCK was 
developed. (See Appendix A for a complete discussion of 
SUBROUTINE PLUCK.) Basically, SUBROUTINE PLUCK uses the re- 
sults of the static program, but causes a crystal to be 
printed out that only contains the interstitial and all lat- 
tice atoms from two planes below the interstitial to the 
surface plane of the crystal and all atoms in vertical planes 
which are within two (or three) lattice units of the vertical 
plane containing the lattice site shared by the split inter- 
stitial. (See Figure 2.) The savings in computer time re- 
sultingfrom the use of this smaller crystal is demonstrated 
by the following comparison. A thirty timestep dynamic run 
with the entire 10 x 10 x 10 crystal (250 atoms) used approxi- 
mately three and one half minutes of computer time. A thirty 
timestep dynamic run using ae ae SX erystal (60 atoms) 
used slightly less than one minute of computer time - a 71% 
savings in computer time: 

Most of the dynamic runs made'during the course of 
these investigations utilized the 7 x 5 x 7 "PLUCK" crystal. 


Final confirmation of minimum energies was determined using 


the entire crystal. 
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2. The Potential Function 

As mentioned previously, the many-body interaction 
which characterizes actual lattice dynamics is approximated 
in computers by many two-body interactions. These two-body 
interactions are represented in the computer by some type of 
central, pairwise potential function. Various types and 
combinations of potential functions have been considered for 
use in computer simulation investigations of lattice dynamics. 

The choice of the potential function must be made 
with due consideration to the range of applicability of the 
potential function, the correlation of potential function 
parameters with observable properties of the material being 
investigated, and the amount of computer time required for 
calculations using that potential function. The potential 
function used in these investigations was the composite 
Born-Mayer potential and Girifalco and Weizer Morse potential 
used by Tankovich [3,16]. This composite potential is con- 
structed as follows: 

a. Region 1-(г < 1.58) 

The atom-atom interaction at close approach is 

represented by a Born Mayer repulsive potential of the form, 


СЕ = exp (A+B rij? (1) 


where vas is the interaction energy between particles i and 
Hunc Eqs is the distance between particles i and j. 
b. Region ЗЕ ЛОВА =< Ir = 2.028) 
This portion of the potential fünction-is obs 


tained by computing the best cubic equation between the value 
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of the Born Mayer potential at 1.5% and the value of the 
Morse Potential at 2.0%. 


C. Region Gon бот = 5.382) 


For equilibrium and greater separations, a 
Morse potential of the form, 

e = IS 2 P (2) 
where IE is the interaction energy between particles i and 
j, D is the dissociation energy of the particles i and j, rij 
is the distance between particles i and j, and Lo is the 
equilibrium separation, is used. The Morse potential was 
computed so that the tail of the function was truncated to 
zero at 5.38%. This effectively meant that atoms out to the 
fourth nearest neighbor were included in calculations of 
interaction energies. Girifalco and Weizer [7] had pre- 
viously determined Monee potential parameters for tungsten 
and other elements which included interactions out to the 
150th nearest neighbor. Use of the complete function, 
however, would have required an excessive amount of computer 
Lime for Calculation. Additionally, contributions to the 
interaction energy of all atoms beyond the fourth nearest 
neighbor is essentially insignificant for our calculations. 
Many computer simulations of lattice dynamics of body- 
centered cubic materials have utilized potential functions 
which only included first and second nearest neighbors with 
satisfactory results [6,11]. To check the adequacy of this 


potential function in describing this crystal system, the 
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largest binding energy observed in the crystal model 
(-8.283 eV) may be compared with the experimentally de- 
termined heat of sublimation for tungsten (-8.8 eV) reported 
by Harrison and Magnuson [17]. This agreement within 5.9 
percent was considered satisfactory for these simulations. 
3. The Timestep 
a. General Discussion of the Timestep 

Most of the early simulations of lattice dynamics 
used a central difference method as the numerical procedure 
of integrating the equations of motion of the atoms in the 
crystal. See Gibson, et.al [4] or Gay, Effron and Harrison[10] 
for an explanation of the central difference method of numer- 
ical integration. The investigations reported in this paper, 
however, used the average force method which is completely 
described by Gay, Effron, and Harrison [10]. Inherent to 
both numerical methods of integration is the replacement of 
time derivatives in the equations of motion with a finite 
time difference, i.e., the timestep interval. As mentioned 
previously, the strongest interatomic interaction places the 
greatest demand on the system; consequently, the timestep 
duration is usually determined through consideration of the 
strongest interatomic interaction of the system. These in- 
vestida tios followed the procedure of Gay, Effron, and 
Harrison [10] and determined the timestep duration by a con- 
sideration of the maximum displacement that the most ener- 
getic atom of the crystal was allowed to move. This parameter, 
referred to hereafter and in all computer programs as DTI, is 


determined as follows: 
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(1) The equation of motion of Де atom of the 
crystal during timestep interval AT can be written in the 


förm 


x. (t*AT)- x; (t)*[v; (6) + (F, )AT/2m] AT (3) 
or in the equivalent form 


i (v. (F; )AT/2m) AT. (4) 


(2) Rearranging terms of equation (4) yields 
nue Ax, / (v ,+ (F, )AT/2m) (5) 


where AT is the timestep interval, Ах. 1s the displacement 
or the 45% atom during the timestep, V; 15 the velocity of 
МЕП ach 
the ı atom, and (F, ) is the average force on the 1 atom. 
From equation (5) it can be seen that the 
timestep interval is a function of both the kinetic energy 
nac he force. Rather than solving this quadratic equation 
for the timestep interval, the average force method considers 
separately the cases when energies dominate forces 
(v; 22(F; )AT/2m) and when forces dominate energies 
( (F; »4T/2m >> vi). Solutions for each of these cases yields 
a different value for the (next) timestep interval. 
In energy dominant cases, 
or, expressed in the form of energies, 
2T ) ^4 (6a) 
AT = AX; (m/ me a 
and, in force dominant cases, 


АТ = (2m Ax,/ (F,)) 5, (7) 
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where Т of equation (6A) represents the kinetic energy 

of the atom of the crystal with the greatest kinetic energy, 
and (Fi) of equation (7) represents the average force on 
the atom with the maximum force, the AX. of each equation 
becomes DTI. 

In these investigations, DTI was set at the 
beginning of the program. Ideally, the proper choice of 
timestep duration for the first timestep and the proper 
choice of DTI would lead to a smooth movement of the inter- 
stitial to its equilibrium position in a manner similar to 
the movement of a critically damped oscillator. 

It was anticipated that energies would 
dominate forces in early timesteps which would lead to 
timestep determination from the energy equation, equation 
(6A). At some point in the crystal relaxation procedure, 
energies should have been dissipated to the point where 
further timestep determination would become force dependent 
(equation (7)). 

b. The Average Force Method and Timestep Determination 

In the average force method of integration of the 
equations of motion, velocities (i.e., energies) of, and 
forces among, all atoms of the crystal are computed with the 
atoms in their initial positions. Based on these forces and 
energies and the initial timestep duration, new positions 
for all atoms of the crystal are computed. The forces at the 
new positions are then averaged with the forces at the ori- 
ginal positions to determine the average force, and hence the 


final positions of all atoms at the end of the first timestep. 
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In the meantime, the maximum force exerted on an 
Sensor tne crystal in its original position and the maximum 
s ex rto n an alom of the crystal in its final (aver- 
aged) position are used with the DTI and equation (7) to 
Calculate two possible alternatives for the next timestep 
interval. Likewise, the energies of the most energetic atoms 
in both original and final positions are used with DTI and 
equation (6A) to calculate two other possible alternatives 
for the next timestep interval. These four alternatives are 
then compared, and the smallest is chosen as the next time- 
step interval. 

c. Procedures Used to Determine DTI 

Nine l IRurtilized a constant DIL in. all of his 
computations. Since the movement of an interstitial to an 
equilibrium position cannot, in general, be characterized by 
a small range of energies and forces and since DTI should be 
closely correlated with energies and forces at least in ap- 
pmepriate regions, the use of a constant DTI in all calcu- 
lations made the initial choice of DTI extremely critical. 
Success could only be attained by resorting to small DTI's 
and concomitant excessive program run times (> 100 timesteps). 

Tankovich [3] obtained more satisfactory results 
b successively decrementing DTI for each timestep. This 
procedure allowed long timesteps in early portions of a run 
when the interstitial was far from its equilibrium position. 
As equilibrium was approached, the decrementing process had 
progressed to the point such that significantly smaller and 


smaller timesteps were taken allowing a smooth arrival at the 
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equilibrium position. Additionally, this was accomplished 
at a considerable savings of computer time (~ 30 timesteps). 

Although this decrementing process for DTI pro- 
vided considerable improvement, on occasion, the final few 
emesteps used such a small DTI that practically no atom 
movement was discernible. During these investigations this 
Situation was alleviated by incorporating a minimum DTI into 
the decrementing process. This insured that atom movement 
was still discernible near equilibrium and assisted in 
guarding against possible false assumption of equilibrium 
because of the relatively small movement observed under the 
continuous decrementing process. 

In an attempt to more fully understand the 
mechanisms of the static solution, the computer program was 
adjusted so that the maximum force, the atom upon which this 
force was exerted, and the four "new" possible timesteps for 
each timestep interval were printed out after each run. It 


Was observed that the timestep calculations based upon "new" 
and “old" energies were overwhelmingly the basis for timestep 
determination. To insure that the force dependence of the 
timestep was not sele unduly disregarded, the program was 
adjusted so that DTI was determined solely as a function of’ 
the minimum of the two forces. These "force calculations" 


of equilibrium positions agreed with "energy calculations" 


within 0.03 lattice unit. 
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B. THE PROGRAMS 
Although the basic computational procedures contained 
in the computer programs for the static and dynamic portions 
of the problem were essentially the same, an understanding 
of the subtle differences between the programs and an ap- 
preciation for several computational tools and procedures 
used in the programs have to be gained prior to further dis- 
cussion of the actual investigation and results. 
IS The Static Program 
In the static program, a tungsten crystal of 
appropriate size was created in the computer. An inter- 
stitial atom was then implanted at a chosen site within 
the crystal. Potential energies and mutual forces of all 
atoms were computed. The crystal was then allowed to relax 
in appropriately chosen timesteps. At the end of each time- 
step an energy dissipation mechanism was introduced in the 
form of a predetermined damping factor which was used to 
decrease all velocity components of the atoms of the crystal. 
The next timestep interval was then computed, and the process 
was repeated until an indicated number of timesteps had been 
completed. If an equilibrium position had been reached, 
positions and energies of all atoms in the crystal, including 
the interstitial, and other pertinent data could be punched 
out on cards for later use in the dynamic program. 
a. Equilibrium Positions of Interstitials 

In interpreting the results of the static simu- 

lations it was necessary to arrive at some criterion to use 


as a determining factor for the interstitial's final 
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equilibrium position. It was expected that the crystal 
wonllafrelax around the split interstitial site forming a 
"pocket" within the crystal which would interrupt the 

e eeit oE the lattice. Stability in this configuration 
was determined to have been reached when the atoms of the 
deformed crystal, including the interstitial, had been 
allowed to alles: (i.e., adjust to the presence of the 
interstitial) to the point where their kinetic energies 

were all below thermal.  (« 0.025 eV) 

If more than one possible equilibrium position 
met this criterion, it was felt that different positions 
within the "potential well" of the equilibrium site were 
probably being observed. The position in which the inter- 
stitial atom had the smallest amount of potential energy was 
then chosen as the equilibrium position for the lattice site 
under investigation. 

Den nalng Of Oscillations Near Equilibrium 

While performing the "force calculations", some 
runs began with the interstitial moving toward an apparent 
equilibrium until an oscillation, or "rattle", developed 


about the suspected equilibrium position in the "z" direction. 
It was first confirmed that this "rattle" was solely in the 
Pee ctrection (amen, no significant “x! or "y" displacement) 
by extending the computation from 30 to 200 timesteps. It 
was shown that in a (100) direction in a body-centered cubic 
a narrow potential energy minimum was observed at (100) 


planes. The region of the (110) line between the octahedral 


void and the reference lattice site was contained within 
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this potential energy minimum. With this information as 
ISE catión, TE was determined that the velocity of the 
intérstitial in the "z" direction could be completely damped 
at the end of each timestep to hasten the determination of 
the true equilibrium position. This technique was useful 
in restricting computer run time whenever an obvious "rattle" 
in the "z" direction developed. 

2. The Dynamic Program 

The initial step of the dynamic program was the re- 
creation in the computer of the tungsten crystal, including 
the interstitial, after relaxation of the crystal had taken 
place and the interstitial had come to an equilibrium posi- 
tion. This was accomplished by utilizing the output of the 
Static program as the input to the dynamic program. 

Before continuing with the problem, it was necessary 
to realize that the minimum energy required for the inter- 
stitial to escape from the crystal would be a function of 
the path of escape. To account for this "direction depend- 
ence", impact points were chosen in the surface plane of the 
crystal perpendicular to the direction of escape. The nega- 
tive "y" direction was always used as the direction of 
escape. Energy was then imparted to the interstitial which 
was, in ewEn; "armed" atwa specific impact point. The inter- 
stitial was then allowed to travel through the crystal using 
a timestep procedure based on a constant DTI. (The timestep 
Procedure in the dynamic program is simpler than in the 
static program since energies are dominant throughout.) Each 


impact point (i.e., direction of escape) was subsequently 
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tested in the same manner using the same initial energy. 
If the interstitial escaped from the crystal in any di- 
rection tested, it was assumed that the initial energy 
applied to the interstitial was greater than the binding 
energy for that particular ion in that particular equili- 
brium position. In this case, the initial energy was de- 
creased and another survey of the impact points was taken. 
When the minimum initial energy which still allowed the 
interstitial to escape was determined, the binding energy 
e that particular ion in that particular equilibrium 
position was known. 
3. Equilibrium Sites Viewed As Potential Wells 

One means of modeling equilibrium positions of 
interstitials in a crystal lattice is to consider each pos- 
sible equilibrium position as a three dimensional potential 
well. A foreign atom migrating through the lattice which 
reaches the confines of one of these potential wells with 
sufficiently small energy would fall into the potential well 
and remain there. In conjunction with this research, the 
character of these potential wells was also investigated. 

Two different approaches were taken to investigate 
the potential well aspect of the problem. First, the static 
program was used to investigate behavior of interstitials 
implanted at various positions around a previously deter- 
mined equilibrium position in a perfect lattice to determine 
which of these positions sought equilibrium. The second 
method offset the interstitial from its equilibrium position 


in the_relaxed crystal and then allowed the crystal to undergo 
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further relaxation to determine whether the offset inter- 
stitial would seek the original equilibrium position. 
4. Determination Of Possible Interstitial Sites 

In a body-centered cubic material there are twelve 
possible locations for a (1105 split interstitial. In an 
infinite crystal these sites are equivalent; that is, no one 
site can be distinguished from another. When a finite crystal 
is considered, the presence of a crystal surface allows 
identification of three distinct interstitial sites. An "A" 
Site is located in a (110) plane which contains the lattice 
Site about which the split occurs and is closer to the 
crystal surface than the shared lattice site. A "C" site 
is also in a (110) plane containing the shared lattice site, 
but is located below the shared lattice site. A "B" site 
is located in the same (100) plane parallel to the surface 
that contains the shared lattice site. 

a. Implantation Procedure 

In the static program the interstitial was 

initially implanted in an offset position in the direction 
Of the expected equilibrium position. This procedure al- 
lowed a smooth movement toward the suspected equilibrium 
position in a minimum number of timesteps. For the heavier 
interstitial atoms (xenon and tungsten) the lattice atom 
was also offset in its direction of suspected movement as 


a result of the presence of the large interstitial. 





IV. PRESENTATION OF DATA 


A. STATIC SIMULATION 
M initial Simulations 

After the minimum DTI procedure had been incorporated 
into the DTI decrementing process, investigations of tungsten 
interstitials in a tungsten crystal were made. In particular, 
a comparison between tungsten movement under the influence 
of an attractive potential and tungsten movement under the 
influence of a repulsive potential was sought. Equilibrium 
positions determined separately with these two potentials 
gave agreement within 0.02 lattice unit. 

Investigations were then extended to verify equili- 
brium positions near lattice sites 89 and 64 which had pre- 
viously been obtained by Tankovich [3]. Positions for argon, 
neon, krypton, and xenon were examined for each site, and 
results agreed with Tankovich's data within 0.02 lattice 
unit, see Table I. As the mass of the interstitial atom 
decreased, its equilibrium position moved from the shared 
lattice site Tocation previously reported [7,8] toward the 
center of the octahedral void. 

2. Site 39 

Investigations were then directed toward determining 
equilibrium positions of a split interstitial near lattice 
ce O аен 5 one layer below the surface of the crystal. 
Calculations were made separately for argon, neon, krypton, 


and xenon interstitials. At the end of thirty timesteps, 
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TABLE I 


Interstitial and Lattice Atom Displacements from Reference Site 








Interstitial Meerstitial Lattice Ato 


. Site 
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all atoms of the crystal in each case had kinetic energies 
below thermal, which indicated that an equilibrium had been 
reached. An observation was made, however, that tended to 
abrogate this determination. The lighter interstitials, 
neon and argon, exhibited a definite affinity for the sur- 
face of the crystal. In all three possible sites (A,B, and 
C), argon and neon interstitials were found to seek posi- 
tions Significantly (~ 0.1 lattice unit) closer to the sur- 
face than the expected (110) split. Final positions for 
argon and neon interstitials in site 39A were located less 
than three tenths of a lattice unit below the surface of 
the crystal. Similar behavior, but to a lesser degree, 

was also observed with the heavier interstitials, krypton 
and xenon. Additionally, although the kinetic energy cri- 
teria indicated that an equilibrium had been reached, the 
small velocities that were available at the end of each 
timestep were in such a direction as to allow escape from 
the crystal if these velocities could be maintained. In 
fact, although velocities were halved at the end of each 
timestep, a comparison of velocities over the last five 
timesteps of the computer run showed that the negative "y" 
(direction of escape) velocity component of the interstitial 
pemene send of timestep thirty was actually only twenty per- 
cent lower than the same velocity component at the end of 
timestep twenty five. In other words, even with fifty per- 
cent damping applied during each timestep, velocities actu- 
ally decreased by only twenty percent over 5 timesteps. It 


was also observed that during these last five timesteps the 
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maximum movement allowed by the most active atom of the 
crystal (the DTI) esa reached its minimum value, 0.0005 
lattice unit. These observations suggested that the small 
movement allowed during these timesteps combined with the 
BLULOmep factor might be "forcing" the interstitial to ex- 
Bubit equilibrium criteria. It was concluded that equili- 
brium positions near site 39 were unstable, if they exist 
cr all. 
3. Site 14C 

Until these investigations little thought was given 
Ll possibility of an equilibrium position in site 14С. 
With no precise knowledge TEE concerning the actual 
quantitative value of the damping experienced by a foreign 
interstitial implanted in a crystal lattice and with the 
possibility of equilibrium sites near lattice site 39, 
some credence had to be given to the possibility of equili- 
brium positions near surface atoms of the crystal. It was, 
therefore, decided to investigate the possibility of an 
equilibrium position in site 14C. It was postulated that 
in the case of light atoms (neon, for example) the inter- 
L iS] position for the 14€ site should be deeper in the 
crystal than the 39A site, as a result of the greater re- 
pulsion of the lattice atom with which the site was shared. 
Simulations showed that the 14C interstitial site was located 
at a distance of 0.43 lattice unit below the surface of the 
crystal while the 39A interstitial site was located at a 
distance of 0.26 lattice unit below the crystal surface. 


When the displacement of this interstitial site was compared 
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with the 64C interstitial site (0.84 lattice unit below 
the next higher lattice plane), the effect of the surface 
of the crystal on interstitials in close proximity to the 
surface was clearly exhibited. Velocity characteristics 
similar to those of site 39A were also observed ааа лла 
that the amount of damping applied and the DTI could be 
forcing the interstitial to exhibit equilibrium properties 
gm this site. In general, it can be said that interstitial 
equilibrium sites in the first two layers of the crystal 
weil) defined in position if they exist at all. 
4. Force Calculations 

Poemenetoned previously, in the course of these 
investigations it was determined that energies nearly al- 
ways dominated forces for the timestep ranges used in the 
Calculations, and, consequently, the timestep was nearly 
always chosen as a function of the energy of the most ener- 
getic atom. To ascertain whether the use of timesteps deter- 
mined by maximum forces would yield better or significantly 
different results, the static program was modified so that 
the timesteps were determined strictly as a function of the 
maximum force. This was accomplished by first printing out 
the maximum force observed during each timestep. By sur- 
veying these maximum forces, several values of force were 
chosen to be used as test forces. As the program was sub- 
sequently run, the maximum force in each timestep was compared 
with these test forces, and a DTI for that timestep was as- 
signed based on the results of that comparison. This 


assigned value of DTI was then used to compute the next 
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timestep interval. These calculations of eguilibrium 
positions based on forces agreed with previous energy cal- 
culations of equilibrium positions within 0.03 lattice 
mmt. 

It should be realized here that the terminology 
"force calculations" and "energy calculations" do not imply 
any significant difference in method for determining equili- 
brium positions. They are merely two different procedures 
for determining the maximum displacement which will be al- 
lowed Бру ‚Һе most energetic atom of the crystal during the 
next timestep. The importance of this parameter (DTI) is 
in the effect it has in insuring the smooth movement of the 
interstitial to its equilibrium position. 

5. Investigation of Potential Wells 

Viewing interstitial equilibrium positions as po- 
tential wells of varying depths is a convenient means of 
modeling the entrapment of foreign atoms by lattice struc- 
es. In order to minimize the effects of the crystal 
surface on investigations, interstitial site 89C, which had 
exhibited good stability and almost perfect (110) splitting 
in previous testing, was chosen for investigation of the 
characteristics of interstitial potential wells. 

а. Potential Well Studies in a Perfect Lattice 

The first approach to the study of the inter- 

stitial potential well utilized essentially the same method 
that had been originally employed to locate equilibrium posi- 
tions. It was postulated that any interstitial atom that was 


implanted in a perfect lattice at or near the coordinates of 
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the equilibrium position which had been previously deter- 
mined for interstitial site 89C would seek the same equili- 
brium site. Interstitial atoms could be implanted further 
and further from the previously determined equilibrium 
position until they no longer returned to it. In this man- 
ner a map of the size of the potential well around site 89C 
could be obtained. 

The results obtained in investigations of xenon 
in site 89C are presented in Figure 3. In this figure each 
coordinate intersection represents an implantation site 
which was tested. The arrows drawn from these implantation 
Sites indicate the direction of movement of the interstitial 
meen that Specific implantation site. The tip of the arrow 
represents the interstitial position after thirty timesteps. 

In considering the data obtained in these in- 
vestigations, two significant observations were made. First, 
an interstitial implanted at the previously determined equili- 
brium position (coordinates (4.52, 3.48) in Figure 3) moved 
from this position to another position which also met the 
Seiteria for equilibrium. Secondly, all other implantation 
Sites also moved to positions which met equilibrium criteria. 
n analysis of the program print out data provided informa- 
tion of secondary importance. All implanted atoms exhibited 
an initial movement in the general (110) direction even 
though final positions were not necessarily in that direction. 
Similar investigations were conducted using argon and neon 


Macerstitials with similar results. 
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In light of this unexpected behavior in the 
Bar erslarereeschesstability of aforeign atom, inserted 
into the crystal as a replacement for the lattice atom in 
site 89, was examined. Investigations were made using neon, 
argon, krypton, and xenon as the replacement atoms. The 


data from these runs are tabulated below. 


TABLE TI 


Results of Replacement Runs in Site 89 


Replacement (Y Displacement 


Atom After 30 Timesteps Energy (eV) Energy (eV) 
Neon -0.049 0.0000 1.0444 
Argon -0.036 0.0008 SEDIS 
Krypton -0.029 0.0001 325120 
Xenon -0.020 0.0000 1:3297 


(L.U.) {Final Kinetic{Final Potential 


The y displacement values are sufficiently small 
phat site 89 must be presumed to be a stable replacement site. 
This conclusion is further confirmed by observing that the 
potential energies are also significantly lower than those 
obtained for the same atoms when acting as interstitials, 


See Table III. 


TABLE III 


Interstitial Potential Energies for Site 89C 


Interstitial Final Potential 
Atom Energy (eV) 

Neon 4.7 

Argon 10:53 

Krypton MO 

Xenon S 
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Nostüurtlerscheck the. stability of this site, 
a kryptron replacement atom was then initially offset from 
Bue site 89 coordinates, 0.7 lattice unit in the "x" and 
"y" directions, and the program was run again. The krypton 
atom moved precisely along the (110) line back toward site 
89. At the end of thiry timesteps, the krypton atom was 
Шосавей O.4 lattice unit in the "x" and "y" directions, 
"x" and "y" velocity components were still in the direction 
of site 89, and potential energy had decreased from 54 eV 
in the initial offset position to 8.7 eV at the end of 
thirty timesteps. 
(1) Interpretation of Results 

Interpretation of the results reported in 
the previous section led to a re-examination of the concept 
of the equilibrium positions of interstitials and raised a 
question concerning the validity of solely using relaxed 
@rystals for determinations of equilibria. 

The equilibrium positions obtained by using 
a perfect crystal were a function of the implantation site. 
This dependency of the final equilibrium position on the 
implantation site suggested that the movement of the inter- 
stitial to an equilibrium site can not be explained in terms 
ОЕШ Шс ес сзеншесите “forcing the interstitial to its 
lowest energy configuration. Rather, these results sug- 
gested that the actual mechanism is a process of local 
"liquefaction" of the lattice as it adjusts to the presence 
of the interstitial combined with interstitial movement caused 


Il wEnteractron with the lattice atoms. This, in turn, 
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suggested that a more accurate means ош оса елла еамі1 
brium positions might be to utilize a crystal which had 
already been allowed to "relax", or adjust to the presence 
ОЕК Ле interstitial. (See next section. ) 

The movement of atoms from nine different 
implantation sites in an area two tenths of a lattice unit 
per side to nine different positions which exhibit equili- 
brium criteria indicated that equilibrium positions are not 
precise positions coordinate - wise within the time range 
of these computations. 

b. Potential Well Studies in a Relaxed Lattice 

The conclusions reached as a result of studies 

of potential wells in a perfect crystal prompted similar 
studies in a crystal which had been allowed to adjust to 
the presence of the interstitial. An argon interstitial 
emo interstitial site 89C were chosen for investigation. 
The relaxed lattice was obtained by using the results of the 
static simulation of argon in site 89C which had first in- 
dicated that an equilibrium position had been reached. These 
results were read into the computer as initial positions for 
the lattice atoms E Bie wimeerstitial. "The interstitial 
was then offset from its position and a new static simu- 
lation was performed. An array of sites was chosen for in- 
vestigation in this manner. The results of these simulations 
with various interstitial offsets are shown in Figure 4. 
The representation in Figure 4 is analogous to the repre- 
sentación in Figure os. fhe start point for each run is 


numbered (1-9 and A,B). 
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The interstitial movement depicted in Figure 4 
was encouraging in many respects. An argon interstitial 
placed successively at each initial site (numbered 1-9 in 
Figure 4) appeared to head for roughly the same region of 
Mme crystal. The potential energy of the interstitial at 
the conclusion of each run ranged in value from 8.8 to 9.8 
eV. A comparison of these energies with the energy of the 
interstitial in the initial equilibrium position found in 
the perfect lattice (-16 eV) indicated that this equilibrium 
region is much more stable than the position determined pre- 
viously. While not precisely defined in position, this 
equilibrium region was located in roughly the (110) di- 
rection from the shared lattice site. 

Since this equilibrium region appeared to be 
located nearer to the shared lattice site [аре 89) than 
offset start points 1-9, two additional computer runs (labeled 
A and B in Figure 4) were made with start points bordering 
this equilibrium region. These runs indicated that, even 
in a relaxed crystal, the equilibrium reached was still 
Somewhat a function of the initial position of the inter- 
stitial. In both of these runs, equilibrium was reached by 
movement in the (110) direction, and the potential energy of 
the interstitial after thirty timesteps (~ 10 eV) was lower 
than the potential energy originally determined in the per- 
fect lattice. 

(1) Interpretation of Results 

The results obtained in ine relaxed lattice 


seemed to reinforce the conclusions drawn during the study 
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of potential wells in the perfect lattice. Basically, the 
relaxed crystal studies gave a picture of equilibrium posi- 
tions much more in tune with what might reasonably be ex- 
pected in nature. It seems clear that the equilibrium site 
of an interstitial is an ill-defined region in the general 
Ш) direction from the shared lattice site. Determination 
of this equilibrium region cannot be made based solely on 
the kinetic energy of the interstitial; potential energies 
must also be considered. 

Р This equilibrium region could have been 
postulated from a consideration of the rms displacement of 
a tungsten atom in a tungsten lattice. Houska [18] has 
measured this rms displacement to be 0.049 lattice unit at 
300°K. The ае Еее activity of the lattice atoms, 
then, can be expected to cause the equilibrium positions of 
interstitials to fluctuate over some region. This equili- 
brium region might best be described as roughly a cylinder 
whose axis 1S in the (110) direction and whose height and 
radius are some function of the relationship between inter- 
Aleta mass, Lattice atom mass, and the interatomic poten- 
tial function., The equilibrium positions determined in these 
investigations appeared to be centered in the (110) direction 
at the intersection of the (110) and (100) planes through 
the lattice atom and covered a section of the (110) line 
on the order of 0.2 of a lattice unit long. This seems to 
be a reasonable range of equilibrium positions for the re- 


latively light argon interstitial in a tungsten crystal. 
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The results of these investigations have 


indicated the usefulness of conducting implantation studies 
in relaxed crystals. The movement of interstitials to equi- 
librium in the relaxed crystal was much less dependent upon 
initial positions, and all sites investigated showed a pre- 
ference for positions in the (110) direction. Such was not 


the case in perfect lattice studies. 


B. DYNAMIC SIMULATIONS 

Concurrent with the investigations into the potential 
wells of the equilibrium sites, dynamic simulations were 
conducted using the argon interstitial in site 64A. Since 
Tankovich [3] had investigated an argon interstitial in 
Site 89C and had determined that the binding energy of this 
Site was in the range 2-4 eV, 4 eV was chosen as the initial 
energy for the dynamic tests. It was determined that the 
binding energy of the argon interstitial in site 64A was 
between 0.02 and 0.04 eV. At 0.04 eV eScape of the inter- 
stitial was noted for several of the twelve impact points 
tested. At 0.02 eV the interstitial first moved slightly 
toward the surface of the crystal, then changed directions 
and moved to a position deep inside the crystal while gaining 
considerable energy. This phenomena had been seen and in- 
terpreted previously as the expected behavior of a stable 
Site in which the interstitial has been required to move 
too large a distance in one timestep. Since the timestep is 
constant (0.1 lattice unit) in the dynamic simulation, an 


interstitial that is oscillating in a stable configuration 


43 





Pale required to move out of its stable position. When 
this movement is directed into the crystal, it can be assumed 
that the original position was stable. 

This binding energy falls in the range of the first 
desorption peak measured by Kornelsen and Sinha (see Figure 
Re. 1). This result indicated that interstitial sites 
formed with lattice atoms in the third crystal plane (cor- 
responding to the y = 2 plane of these calculations), such 
that the interstitial site was in the (110) direction above 
the lattice atom, are the sites nearest the surface of the 
Crystal stable enough to entrap measurable amounts of argon. 
This excellent concurrence with experiment further sub- 
stantiated the conclusion drawn earlier in these investi- 
gations about the doubtful character of positions near sites 


39 and 14 which had exhibited equilibrium criteria. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


During the course of these investigations the possibility 
that the DTI decrementing procedure and the procedure for 
incorporating damping into the problem might be forcing in- 
terstitials to exhibit equilibrium properties arose. At 
first, the DTI was thought to possibly be restricting atom 
movement to such an extent during the last few timesteps of 
a calculation (i.e., when DTI had reached its minimum value 
-0.0005 lattice unit) that atom movement and, consequently, 
kinetic energy could no longer be used as equilibrium cri- 
teria. This became more significant when the small DTI and 
the damping factor were considered together. The damping 
L etror seemed particularly appropriate for scrutiny near the 
A tace ot the crystal and when DTI was small, since the pro- 
bability of collisions and close approaches and the subsequent 
damping of atom motion could be expected to decrease in both 
instances. The excellent agreement with experiment of the 
results of the dynamic and static simulations suggests, how- 
ever, that the DTI and damping procedures employed might be 
adequate for these simulations except when atom movement is 
Ee uurbxeboscocto the surface of the crystal. It is 
recommended that future investigations explore the use of a 
damping factor which varies with decreasing DTI and movement 


toward the surface of the crystal. 
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These investigations have shown the value of using a 
relaxed crystal in the computer simulation of lattice dynamics. 
The relaxed crystal simulation indicated that: 

(1) The equilibrium position initially sought by an 
шибскселЕта! is a function of Implantation position. 

(2) The mechanism associated with the establishment of 
an equilibrium position seems to be a combination of inter- 
atomic interaction and local liquefaction of the crystal 
са иге in the vicinity of the interstitial. 

(3) The character of equilibrium potential wells is more 
readily observable in a relaxed crystal. 

(4) The actual equilibrium position of an interstitial 
Seems to be some region in the general (110) direction from 
epe reference (shared) lattice site. 

Again, however, for determining binding energies, the 
original procedure for determining equilibrium positions may 
yield adequate results. Such was the case for the dynamic 
Simulations reported here. It is recommended that replace- 
ment interstitials and other split interstitial sites undergo 
dynamic simulation in an attempt to correlate other desorption 
peaks as determined by Kornelsen and Sinha [1] with inter- 


stitial sites. 
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APPENDIX A 
SUBROUTINE PLUCK 

Subroutine Pluck was developed to allow the use of a 
crystal smaller than the original model in dynamic testing. 
This subroutine was developed so that the only pieces of 
information required as input to the subroutine were the 
number of crystal planes desired in the "x" and "z" 
directions and the number of the lattice site desired at 
the center of the PLUCK crystal. Once the size of the 
PLUCK crystal is decided upon, the subroutine stores the 
original numbers of the atoms in the PLUCK crystal in an 
way. This allows reference to any atom by its original 
number throughout computation. The atoms of the PLUCK 
crystal were numbered consecutively, and the number of atoms 
in the PLUCK crystal was assigned a variable name. In this 
manner, a minimum of adjustment was required in the dynamic 
program when the PLUCK crystal was used. 

SUBROUTINE PLUCK is included in this appendix in its 
most general form. Parameters are listed below for two 
different sizes of PLUCK crystals. Sizes refer to the 
number of "x" and "z" planes in the PLUCK crystal. "Y" planes 
from the surface of the crystal through the two planes below 
the plane of the lattice site under investigation are always 


included. 
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PLUCK PARAMETERS FOR DIFFERENT SIZED CRYSTALS 


krystal sizes refer to number of planes of atoms in x and 


z directions) 


Parameter 


IXNEW 
IZNEW 
NM1 
NIl 
NX1 
NII31 
NII41 
NIINCI 
T TITINC1 
NMINC1 
NXINC1 
NMINC2 
NX2 
NM2 
N12 
ug 32 
NII42 
NIINC3 
IIINC2 
NMINC3 
NXINC2 
NMINC4 


maz 7 Crystal 
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10 
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ALPHA: . 


BSAVE: 


BIND: 


BMAS: 


BULLET 


APPENDIX B 

COMPUTER PROGRAM GLOSSARY 
Input Morse potential parameter. 
Target Mass/ (target mass + bullet mass); 
distributes potential energy between target and 
bullet. 
Negative of the total potential energy (TPOT) 
at time zero. 
Mass of bullet in AMU. 


Alpha-numeric array for point defect material. 


CFO, CFl, CF2: Force parameters of cubic fit between Morse 


ева, 
CGD1, 


"GEI, 


CGB2: 
CGD2: 


CGF2: 


and Born-Mayer functions. 
Morse potential parameters. 
Morse potential parameters. 


Morse force parameters. 


error CPI, CP2, CP3: Potential parameters of cubic fit 


between Morse and Born-Mayer functions. 
10 


ЕСО СУБ X 10 ©: Converts lattice units to meters. 

EYE: 1.6 X om Converts electron volts to Joules. 

CVED: CVE/CVD: A ratio used to avoid repeated division. 

CYM: 1.672 X Troma Converts atomic mass units to kilograms. 


CVR: LU in angstroms; converts lattice units to angstrom 


units. 


DIX, DIY, DIZ: Displacement coordinates for location of 


DCON : 


interstitial from reference atom, NVAC. 


Input Morse potential parameter. 
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fÉ a 








DOTIPF: 


DDTII: 


DFF: 


DIST: 


The minimum value that DTI is allowed to assume. 
The initial decrement of DTI. 

ROE-DIST. The distance closer than ROE that an 
atom is to the primary. 


Distance between any two atoms. 


DLPE: TLPE-TLPE®@: The change in total local potential energy 


since time zero. 


DEX, DRY, DRZ: X,y,Z components of DIST. 


DT: 


bee, DTEJ: 


DEE, ОТЕ1: 


DTI: 


Length of a timestep in seconds. 

The two possible alternatives of the timestep 
computed from maximum energies. 

The two possible alternatives of the timestep 
computed from maximum forces. 

Number of lattice units most energetic atom may 


move in one timestep. 


Pees) DT2 (I), DT3 (I), DT4 (I): Vector arrays which save 


DTSTEP (I): 


the possible choices of timestep determined in 
the "energy” method. 
Vector array which saves the timestep interval 


chosen for each timestep. 


DTOD: DT/CVD -- A ratio used to avoid repeated division. 


DTOM:DT/PTMAS -- A ratio used to avoid repeated division. 


pau s bDYOUbAXG):"change in position of ith atom from 


EMAX 


EV: 


EVR: 


initial position at time zero. 
The maximum energy encountered in any cycle. 
Primary energy in electron volts. 


Primary energy in kiloelectron volts. 
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EXA, EXB: Input Born-Mayer potential function parameters 
for the target. 

F2: ` Square of the force on a specific atom. 

FA: The component force increment on an atom. 

FDTI: DTI X CVD: A parameter used to determine DT by 
maximum energy method. 

FM: A small number used in checking potential energy 
Zero point. 

FM2: FM squared. 


1 


ЕМАХ Maximum total force on the most stressed atom 


in the crystal. 
FOD: FORCE/DIST: A ratio used to avoid repeated division. 
BORCE: Numerical value of the force function with a 
variable parameter. 
FORMAX (I): Vector array which saves the maximum force in 
each timestep. 
mae EY (I), PZ(I): x,y,Z, components of total force on 
an tacon: 
PAX: Born-Mayer force function parameter. 
HBMAS: 4 BMAS: A ratio used to avoid repeated division. 
HDTOD: % ртор: А ratio used to avoid repeated division. 
HDTOM: % DTOM: A ratio used to avoid repeated division. 
HDTOMB: 5 DTOMB: A ratio used to avoid repeated division. 


HTMAS: 4 TMAS: A ratio used to avoid repeated division. 


Il: Vertable in -eubic frt subroutine. 
13: variable in Cubic fit subroutine. 


IDEEP: First fixed layer. 
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IFMXAM (I): Vector array which saves the atom number which 


experiences the maximum force. 


EHI: Alpha nümeric array for program title. 

IH? : Alpha numeric array for Morse function parameters. 

«BEP : | Alpha numeric array for bullet element. 

IHS: Alpha numeric array for type and orientation of 
crystal: 

DHT: Alpha numeric array for target element. 

ШІ : Number of atoms in а crystal using subroutine 
PLUCK. 

ILAY: Number of free (mobile) layers. 

IN: Odd-even integer used to determine atom site 
establishment. 

ЇР: Subscript value of atom. Void in subroutines 
STEP and ENERGY. 

Tipi: Parameter that determines whether or not a self 
defect is to be given a repulsive potential or 
a composite attractive - repulsive potential. 

SHUT: A parameter used БОК ЫК асбуп che program. 

ш: Unscaled fixed point x coordinate used in 
lattice generation. 

ITT: Odd-even integer used to determine atom site 
establishment. 

JT VPE: Parameter used to determine the type of point 


defect: vacancy, interstitial, or replacement. 
Ша ту то ANDES IO x yo Z planesiof cyrstal. 
IXNEW, IYNEW, IZNEW: Number of x,y, and z planes in the 


PLUCK crystal. 
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J2: 


JU: 


JT: 
JTS: 
JTT: 


KEEP (I) 


MENT (I): 


LD: 


LL: 


LOCAT (K): 


LS: 


MCRO: 


NEW: 


NNUM (1): 


Variable in the cubic fit subroutine. 

Parameter in the BCC(111) lattice generation 
subroutine. 

Unscaled y coordinate used in crystal generation. 
Variable used to establish atom sites. 

Variable used to establish atom sites. 

Vector array which saves the original atom 
numbers of the PLUCK crystal. 

Final K in LOCAT (K) assigned to an atom. 
Unscaled z coordinate used to establish atom 
site. 

Used to identify an ith atom which is not 
included in calculations. 

The highest numbered atom in the mobile layers. 
The highest numbered atom in the entire crystal. 
Dimensioned variable that remembers the numbers 
of the atoms within a radius ROEL of the primary 
atutime zero: 

Variable associated with each of the nine lattice 
generator subroutines. 

One number higher than the order of the fit 
between the Born-Mayer and Morse potentials, 
always 4 in this simulation. 

Data output increment, in numbers of timesteps. 
Parameter used to determine whether or not atom 
numbers have been stored in LOCAT (K). 

Vector array used in Subroutine PLUCK to re- 


number atoms. 
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NPAGE: 


NRUN: 


NS 


NT: 


NTT: 


NVAC: 


PAC: 


PBMAS : 


PEXA, PEXB: 


PIB C з 
РЕХА: 
РКЕ (Т): 
PLANE: 
POT: 


PPE (I): 


BREINT (I): 


PPENCK: 


PPENEG (I): 


PPEPCK: 


PPEPOS (I): 


Page numbering variable. 

Parameter used to determine whether or not to 
read additional data cards. 

Initial print statement timestep number. 
Timestep number. 

Timestep number limit before shutdown. 

An atom number used to establish point defects 
or used as a reference point for interstitial 
placement. 

Parameter for bullet force function correction. 
Primary mass in kilograms. 

Input Born-Mayer potential function parameters 
for the bullet-target interaction. 

Primary force function evaluated at ROE. 
Primary force function parameter. 

[mestre energy or. the ath atom. 

Alpha-numeric array for lattice orientation. 
Potential energy between two atoms. 

atom. 


Potential energy of the ith 


Vector array that saves the difference in 


potential energy before and after implantation. 
Potential energy check value which determines 
potential energy decreases which will be printed. 
Vector array which saves potential energy 
decreases. 

Potential energy check value which determines 


which potential energy increases will be printed. 


Vector array which saves potential energy in- 


creases. 
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PPESAV (I): 


Vector array which saves the initial potential 


energy of lattice atoms before implantation. 


PPKEEP (1): Vector array which saves potential energy dif- 
ferences between perfect crystal and relaxed 
euystal, 

РРТС: Primary potential function evaluated at ROE. 

ENS (I): Total energy of the ith atom (potential + 
Kinetic). 

PTMAS: Target mass in kilograms. 

RE: Input Morse potential parameter. 

RO: Spacing Constant in BCC(110) lattice generation 
subroutine. 

ROE: Nearest neighbor distance. 

ROE 2: ROE squared. 

ROEA: Maximum cut-off for Born-Mayer potential. 

ROEB: Minimum cut-off for Morse potential. 

ROEC:' Maximum cut-off for Morse potential. 

ROEC2: ROEC squared. 

ROEL: Ra MUS inside of which local potential energy 
is found. | 

ROEL2: ROEL squared. 

ROEM: POE DTI Region in which modification of 
repulsive force must be made. 

RT), RY (D RZ(I): SS y,2/ coordinates of an ith atom 


RXI (I), RYK(I), 


at any time. 
ЁК (ПОО Л У,27, Coordinates of an ith atom's 


їн та ро51Е1оһп. 
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ЕС) БК) БАКЕ): x,y,z coordinates of temporary 
Postetonmomean ith atom during force cycle. 

RXNEWI, RYNEWI, RZNEWI: Vector arrays which contain the 
X,y,Z, coordinates of the atoms of the PLUCK 
crystal. 

RXSAVE, RYSAVE, RZSAVE: X,Y,Z, coordinates of the impact 
point in the dynamic progran. 

SAVE: И РОТ. 

БЕ SCY, SCZ: X,Y,Z, coordinate scale factors. 

SSCZ: ПИСЕ ЛЕ ассо сеа for the FCC (111) Lattice 


generator subroutine. 


START: An optional timing variable, not used in this 
simulation. 

SUM: Name lefa recubre fit subroutine. 

TARGET: Alpha-numeric array for target material. 

TSAVE: Bullet mass/ (target mass+bullet mass); 


distributes potential energy between target 
and bullet. 

ТЕ: Total energy of all crystal atoms (kinetic + 
potential). 

TEMP: Temperature of lattice in degrees Kelvin not 
used in this simulation. 

TFAC: A time factor ratio used to determine DT by 


maximum force method. 


ТЕАСВ: TPAC for the bullet. 

THERM: Thermal energy of atom not used in this simu- 
lation. 

TIME: Elapsed problem time in seconds. 
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DOPE: Total local potential energy of atoms within 


a radius ROEL. 


ТІРЕб: ТОРЕ л Cime zero. 

ТМАЅ: Target atom mass іп АМО. 

ТРКЕ: Mane ie energy Of all crystal atoms. 
ШРОТ: Total potential energy of all crystal atoms. 
MISS : Storage variable for velocity components. 


VS (1), VY(1), VZ(1): x,y,z components of ith atoms velocity. 
FV, 2: Unscaled coordinates used in crystal generation. 


XNVAC, YNVAC, ZNVAC: The initial displacement (in LU) of 


atom NVAC. 
YLAX (I): Relaxation in -y direction of ith layer in L.U. 
p losing point form of JTT. 
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APPENDIX C 
ERROR IN THE LITERATURE 
It was discovered during the course of these investi- 
gations that Equation (20) of Ref. 10 (corresponding to 
equation 6A of the report) was incorrect. 


Equation (19) ot Ref. 10 is, 
nu = Ax; /V; 


where AT is the timestep interval, AX; is the displacement 
. th : : : EN 

of the 1 particle, and v; is the velocity of the ı 

particle. 

To express the timestep in terms of the energy of the 
particle with the maximum energy, AX; is defined as the 
displacement of the particle with the greatest energy and 
is replaced by the symbol D. T, is this particle's energy, 


which is the largest kinetic energy at end of each timestep. 


ТГ Т is expressed as 


Eobstrtution into equation (19) yields, 


L 
D/(2 T /m)*, 


AT = 
1 
E un цев со u mu 
vice L 
AT = D (2m/T.) С Equation (20)Ref. 10) 


Tt should be noted that this error has no affect on the 
calculations reported here. D (DTI in these calculations) 


is a parameter which has been specifically chosen in order 


39 





to provide a timestep interval over which there are no 
КИШ Ел Оста torce changes. The factor of 2 which is intro- 
duced in this correction would merely have required a sub- 
sequent alteration of the DTI's chosen so that the timestep 


interval would remain essentially the same. 
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Figure 2. The Function of Subroutine Pluck. 
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64 


Е 





COMPUTER PROGRAM 
FOR STATIC SIMULATIONS 


e 
“шш 
— 
— о 
WO W 
ruil < > 
-LOO > e co 
uc. MGE UI 
KIA tc  €9 €) 
H| 20 UJ с W 
«iex ©) с A Z 
OSM W wea 
4 Z oe - = 
= <0 a 
лула > > ul 
uwuZz>r 0 da + 
> O O < UJ 
E EN E 
Seta eZ 
аа Ge) areal 
OJ = = 4 
Ecd0O0u rp er F 
= — 4 
OrrZ YN nn > 
W=2"Ñ"00 O œŒ 
ru р-а. а ©) 
2MAu 
Cees 2 æ 
LLJU JUL Us = 
gru Un =— J 
WZeOMDY A .J 
SeOcgtl m 0 < 
Ш. æ 2 
Qua bac 3 = U) 
202и = 


Jr QOwu-«U — 
DeilurZrzZ2r> 
Zur Zr < —<r O 
HY III 
@ O> = x Dr. Dr 
QCD ош 2 
C) Oo OL. JI m 
Е ОЕ cv. 
= >Z OU. O GO Oo. 
Teo a 
Бш а О m~ ~A 
non Det N V 
Ia Ш 
eL) EU) 
DAL 
о Ії 
UJ «T co I- 


J IPILUCK FOR USE IN INITIAL DYNAMIC 
оше И ПНЕ  РОТЕМЕТАС WELLS IN THE RELAXED 


Ur OO С) 


CeO 
ш =) Fe 
га по) e y = 


Mee EOS 
R EACH CONFT6= 


LUS 
с 
ИЕК 


XC 
ME 
EN 


M USED E 
NINAS OE 
NAS 


W 


ARRA 
Оа 
ВЕЕШ 


CC OOO 


жж жокк эк Жл ЛО Ух Эк АС >K >K 3 >K Жс Жк Жс кс ж как ок XK 3k 2 К 2 Ж Zx XK 3k XK 2 Xx АК КОЖО КК ЖКЖ X > OK OK OK 


S si ziehe sje sic Жжж Ж ake a ste dete ake ek oe a ae ake ke ke oh te ate a ae sie e o ole ole 


SERRCE CONFIGURATTON! 
xk sesto о ж О ж re К 


es PN Os ASE (S50) ,0TESP(50),0TIVSP(50), 


Z 

Z 

T 
STACC(100),FSSACC(100), FORMAX(100), 


> 2: 

> pr! 

cm mal) 

СЭ СОЕ 

Оса. 

UNUS е 

> > O OO 
> GO OO 
e el AN 
OOo а 6. 
оо UJU 
LALA UN LLik~ 
~ ~ LIJ Y U) 
> >< O. а. 
—-O0.—u-—-O0 oO 

© Оо 

Т гг @ =Z O = Fi = = 
OOO O — O O 
— — ч че — e а 
MMN OY SENnnMw 
ZELEZIUIZAZ* ZE 
LLJUJLLJ Z U PX UJ uu) 
SS a Se x = 
a is AL A LL ач 
OQoOooO0.0-—oOoaoO 


a om 


3k Sk Sk kc ah ae a ae ake ake ake E жож Жуйе Ж жок ж оЖ ate ake ok ok oe otk ac a ok ak ks Жжжж 


Omen CONFIGURATIONS” 
жк жок ai ak a Ee ic ke ae SE aK ae ole zie ck ofc ol oic oleo ol Sc DR AK жж жк ток жжж 


LIONS COO tios 201 
) 


Й 
А 
Т 
STACC(109) )FSSACC(190),FORMAX(100), 


oz 

e ore 

e e” | j) 

ODO 

Ooa 

AN е 

> > O a 
>00 OSO 
e elf el CN 
Q O> oca 
O O< шш 
LALA _ ше 
~ =r L|) HLA) 
> >< Q. H- 
>AULIL-OAO 

O O 

т c) S = = = 
OOOWOROOO 
<-4+-4рь4 че 4 че реа =ч 
(n n Оо) 07) OO) 
== =u >= <= == 


шшш ыш ж шшш 

ccm TIS 

A HA O A LL A — 

ajaja Ma] Jajaja 
acm e 


жж жок ж Ж: OK око he Ae a oos o ccc oc ЗК oc e ooo xe oc sick: oc ole xoc sco oco 5k 2 3k 3k aA a aK OK a 3K 3K >F > ж 


65 





«1 
e > = 
Ly P uL Lu 
Кг e =~ а Ze e 
со e O e >B eN en 
<I О = e @° < LNA O 
pes e > (Y7» > e. JC) UN 
a ~ Fo POn UW OW = N 
< © ae Wy WINE e FONS м, 
> DA ДЕ en e © жо са. ut 
LA ep СС — >< >< LL ec») = 
(С; —<— — = 2. — UY ON ul 
= = OO ONO а ell LL = 
«1 D ~£ ec > e ODAO N 
O mu ж-м ева аі ачи сс 

VY) ы ДЕ OSW KO eee . 
ш e m. OOW UJ etU O 4 == 
er E еШ!) WOO OO. maak e 
co O Az ou ео Сус un 
«1 O Oct m~e O Ok s= «NI 
— N AL NYN XxX ЦУ е СМ е mw 
Om ° — ws CL t >= Ç UJ ON or 
«1 uU») M UY е e Cx. U) ~ мош =z 
=Z uj ac TH meer < LOODU 

Z e meat Om) x *OXNOUOIZZ 
LL ~ re es OOL) Lu — — = .>nN> 

22 O Oe N. > OWWON e e 
WO N nn мә UJ NOOO »-Zz- С) 
Ç -- СУХ v O -uüO.dwcQoule A 
< OÓ > I“ сж“) A >eececlíZin UO 
СС 29) CO mer eQ € e UEONS ey >= 
OW - el) — er CJ *DOUJO Zu o - 
NZ O OJ Ods O OFX + 5 2 Ош 

=! O ND Uy J > NI O € e eCOLLILU LS e 
Z NO-D mm) *. AOIZ ZN% 
Об) —[— e m = 2 Ш -р-цО>жж 
zw X>TM жь-Жж>Ошжа ОО2-+-с23 < 
=> ви CI —«1 QqLe—-PLuadoguouxMI 
Oe Sm en NOS e TRON e nO ALO 
Sie) и СУ-И. СО = + > Ц ОСО — O €N ` — — Z << 

C ZAS еь-.-5 р-5()> 032 5 Ц.3 #3 ез. 
uJuJ QO»?-0O^7FOOQ«U eoOOoooOooo--o 
eoo ORFs om «¿O > DISCO OVVVO 
= М МК МЫМ на Км Мм ММ М UN А 
c vn Z eZee er eS ew*Z<ZzZ ZU.Z ZZ €N Z 
C > QoQOUu--O0--—O»xOuUXDOOoOOO-c 
n < = IS ORNS DOD REST "> 5 За = 
W í =: CONF =з <>. „> > 0X2 ul 
се си Saa IC e] JO» ¿ae Saja (1 aio 
e. < IS et) OL с 

4 чеч — c ocn = 


OOOO 


CER ORNATI STATEMENTS OF ALL READ COMMANDS. 


) 
6.2) 
4X,3F5,312) 


e2 
SF 
+, 


Е 5 
AG 
5 1 
) 


3,F8.4,14,F5.4) 


NITRO 


e N 
SS = 
м am 
- cx = 
+ O 
мч Omm 
а» ы 1 
5 що 
N eD МЈ 
Ф Le e (C) O 
wn en 
a ~ Ce 
= - ad 
< - tU >è 
ME e = => 
= to uu 
C <T 09 
O Oc Ur 
CU) re > 
UJ xz Fm O 
сш ше 
m 24 
ec) — — 
wz. SU 
чє E 
we СО —_ > 
Y) < Od — 
> ka rn 
CONS Oli л 
Б х2 0 
< O eu 
e > < 
LL Z +U i — — e 


ОС) е A eM 
OT NZ 
> FUL UA é =O 
Ar -NZWNUM 
< (CN LULO 
O Mr e< 
0-0 «Ое с 
— rep GQ l е 
YM me + 
c e LL. —J CO) ” 
> uy e <í — >< = 


u IS (P AT STATEMENTS OF ALL WRIT: 
OR 


TT 
OO 00 OO 
mt CJ 


OVO 


pig 


Ud 
Y e 
ON 
mi e 
Ц е 
сп) 
oY) 

e «x 
eL 
Odo 
——0 
иш 
гс, ~ 
e x% 
(yb ex 
eC e 
O N 


Ц cu eme 
MILI = = 


> +m 


t et et Pe |) Se 


° — >< СС 


we 


SUMA 
THAN 


FERRER Oooo 
«t «1 «1 «1 «1 «100 Z C «CD 


== г cm 


— Z s 


cx C. cc CX. ax. xc ac UD UL. «1 cx. vo) 


ro 


-ОО-юсоожосососошхтош 
пиши с“ UL Iu ug u gcua uu cou o 
e 


emer e 


OO NM 
Ct мо г-отоссосо 


DO | QOO QOO OO 
QOO OOO IO O OI O. 


u^ 
о 
JO 
Q^ 


66 





C 
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ТАШУ NAL X INTETAL/EINAL Y IN A NAME РЕ 
PAGHANGE* 7 ) 

Agee a ПИЕ FO. 2: 

СООТ O NN INTTIAL SUMMARY OF POTENTIAL ENERGY CHANG 

КИШ ООО ЛОО АТОМ. РЕ С 

9699 FORMAT (40X,'FORCES AFTER EACH sTEP', //,' TEPIS TOX; 
eee ey tA FMAX 9X, TATOM WITH MAX FORCE',5X, 
IU ORCE Gl)" 7) 

9700 FORMAT (2X,13,6X,7212.4,6X,E12.4,13X,13,12X,E12.4) 

AZ DI] CHECK FOR EACH STEP',//,*'STEP',11X, 
ОЕЕО ҮТ" ОТЕТ ",]]Х,'ОТРЕ1*,11Х,'ОТЕ'*,11Х,'ОТРЕ*,//) 

ERE DOO DI CHECK FOR EACH STEP',//,' STEP' ,11X; 
JEDE age PTF i ts LIX, SDTE*s11XAs*DTFt,//) 

BE AE ATK oe AS BK ae hc SEE NS Жж ot ae К С С С КК К ЖКК ЖЖЖЖ СМ К oe oe I EK Oe a Oe Oe Oi ж кж ж жж 


жжокжож сс Ж жж: УСК 0 К Жжж ж сс ie oe ake NA A A A A A K 


' FORCE CONFIGURATION' 
sk ie I ж жо ж е К жж К Sk SK 3k AK Sk IK ae oe a oe a Sk 3k al 


CHANGE STATEMENT NUMBER 9704 TQ READ 


ENS DN 2X DOT CHECK FOR EACH STEP',//*STEP',11X; 
ШОШ УОИ УКУ Dist ,IIX$,.'DTIV';,//) 


He ak z 3 obs ae SK e E K IK Д АСК к Ж жск: ж He ae i oe ac ac ode a a oe ae ok eK oe SK ae ee e К КУ ж СК ЖК ЭК ЖС 


9705 FORMAT(15,4(3X,E12.4) ,/) 


E 
2 Mie viel Ze eArPPROPRIATE VARTABLES AND VARIABLE ARRAYS. 
START=0.01>*ITIME(XX) 
DO 2 I=1,10J) 
RXK(T1)=0.0 
RYK(1)=0.0 
RZK(1)=0.90 
VX(1)=0.0 
VY(1)=0.0 
М(Т) 59. 5 
PKE(1)=0.0 
РРЕ(1)=0.0 
PTE(I)=9.0 
2 Re2lt1)=0.0 
ISHUT=1 : 
NRUN=0 
C 
С READ INPUT DATA COMMCN TO ALL DESIRED CALCULATIONS. 
READ ( 5,9010) IHl 
КЕДШ 15,9920) IH2YDCON,ALPHA$,.RE,ROEC.ROGEL 
ВЕЦ ОС ОЗО) ВОЕЕЕТ „ВМА „РЕХА „РЕХВ,Т НВ, THERM 
READ I 5,9030) TARGET,TMAS,EXA,EXB, IHT,TEMP 
READ (5,9042) PPEPCK,PPENCK 
C 
Е ИЕ ИЕ ИСО УОЛГА 5 AND DEFINE SCALING FACTORS. 
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